#########################################
## Pratt (2020), Angling for Influence ##
## Replication File #####################
## October 2020 #########################
#########################################

DF <- read.csv("~/Dropbox/DevBanks/Rscripts_2018/Angling2020.csv", row.names=1)

DF_ISQ <- DF[,c("state", "ccode", "year", "institutional_prolif", "prolif_bank", 
                "VotePowerMisalignment", "VotePowerMisalignment_lag",
                "VotePowerMisalignment_otherDBs", "VotePowerMisalignment_otherDBs_lag",
                "Mikesell_Formula_Change", "region",
                "NatIncome_Combined", "Imports_Combined", "Exports_Combined", "Gold_FEX_Combined",
                "time_since_prolif", "time_since_prolif2", "time_since_prolif3",
                "aid_received_aiddata_log",
                "WBvoteshare", "GDP_share", "polity_lag", "aid_given_aiddata_log",
                "aid_given_aiddata_lag_log", "aid_received_aiddata_lag_log",
                "DBs_joined", "GDP_lag_log", "GDPpc_lag_log", "NMC_lag", "tot.trade", 
                "UN_idealpt_USsim_lag", "WB_portsim_st", "portsim_noaid",
                "prolif_join", "prolif_security", "prolif_trade", "regionalDB_voteshare")]

write.csv(DF_ISQ, file="~/Dropbox/DevBanks/Rscripts_2018/Angling2020_data.csv")
DF <- read.csv("~/Dropbox/DevBanks/Rscripts_2018/Angling2020_data.csv", row.names=1)


#################################
## Descriptive Stats & Figures ##
#################################

## Figure 1

v2014 <- DF$WBvoteshare[DF$year==2014]
names(v2014) <- DF$state[DF$year==2014]
v2014s <- sort(v2014, decreasing = T)[1:30]
par(mar=c(5.1, 6.1, 4.1, 2.1))
barplot(sort(v2014s), horiz=T, xaxt="n", yaxt="n", col="lightblue")
axis(side=1, at=c(0,.05,.1,.15))
title(xlab="World Bank Vote Power (% of Total Votes)")
for(i in 1:length(v2014s)){
  par(las=2)
  mtext(side=2, at=seq(0.75, 35.5, length.out=length(v2014s))[i],
        text=names(sort(v2014s))[i], cex=0.7, srt=90)
}
abline(v=.05, lty=3, col="lightgray")
abline(v=.1, lty=3, col="lightgray")
abline(v=.15, lty=3, col="lightgray")



## Figure 2

par(mar=c(5.1, 4.1, 4.1, 2.1))
plot(-1:1, -1:1, type="l", lty=3, xaxt="n", yaxt="n", xlim=c(0,0.14), 
     ylim=c(0,0.14), xlab="Relative GDP", ylab="World Bank Vote Power")
axis(side=1, at=c(0, .02, .04, .06, .08, .1, .12, .14))
axis(side=2, at=c(0, .02, .04, .06, .08, .1, .12, .14))
states <- c("Japan", "China", "Germany", "UK", "France", "India",
            "Italy", "Brazil", "Mexico", "Turkey", "South Africa",
            "Belgium", "Netherlands", "Saudi Arabia")
adj.lab <- c("Italy", "UK", "Netherlands", "South Africa", "Turkey", "Belgium")
for(i in states){
  VS <- DF$WBvoteshare[DF$state==i & DF$year==2014]
  GDP <- DF$GDP_share[DF$state==i & DF$year==2014]
  color <- ifelse(GDP > VS, "red", "black")
  position <- ifelse(GDP > VS, 1, 3)
  points(GDP, VS, col=color, pch=19)
  if(i %in% adj.lab){
    if(i=="UK"){text(GDP, VS, labels=i, pos=4, cex=.85, col=color)}
    if(i=="Italy"){text(GDP, VS, labels=i, pos=4, cex=.85, col=color)}
    if(i=="Netherlands"){text(GDP, VS, labels=i, pos=4, cex=.85, col=color)}
    if(i=="South Africa"){text(GDP, VS, labels=i, pos=1, cex=.85, col=color)}
    if(i=="Turkey"){text(GDP, VS, labels=i, pos=4, cex=.85, col=color)}
    if(i=="Belgium"){text(GDP, VS, labels=i, pos=2, cex=.85, col=color)}
  }
  if(!(i %in% adj.lab)){text(GDP, VS, labels=i, pos=position, cex=.85, col=color)}
}


## When states create alternative institutions, do they have more vote power?
EIBs <- DF$state[which(DF$prolif_bank=="European Investment Bank")]
ADBs <- DF$state[which(DF$prolif_bank=="Asian Development Bank")]
IADBs <- DF$state[which(DF$prolif_bank=="Inter-American Development Bank")]
AfDBs <- DF$state[which(DF$prolif_bank=="African Development Bank")]

mean(c(DF$WBvoteshare[DF$state %in% EIBs & DF$year==2014],
       DF$WBvoteshare[DF$state %in% ADBs & DF$year==2014],
       DF$WBvoteshare[DF$state %in% IADBs & DF$year==2014],
       DF$WBvoteshare[DF$state %in% AfDBs & DF$year==2014]), na.rm=T) # 0.012 avg. vote share in WB

mean(c(DF$regionalDB_voteshare[DF$state %in% EIBs & DF$year==2014],
       DF$regionalDB_voteshare[DF$state %in% ADBs & DF$year==2014],
       DF$regionalDB_voteshare[DF$state %in% IADBs & DF$year==2014],
       DF$regionalDB_voteshare[DF$state %in% AfDBs & DF$year==2014]), na.rm=T) # 0.041 avg. vote share in new IO

## Figure A1
aiddata <- read.csv("~/Dropbox/DevBanks/Data/AidDataCoreFull_ResearchRelease_Level1_v3.1.csv")
aiddata$WB <- as.numeric(grepl("World Bank", aiddata$donor))
yrs <- sort(unique(aiddata$year))
yrs <- yrs[-which(yrs==9999)]
WBcont <- sapply(yrs, function(x){
  sum(aiddata$commitment_amount_usd_cur[aiddata$year==x & aiddata$WB==1], na.rm=T) / 
    sum(aiddata$commitment_amount_usd_constant[aiddata$year==x & aiddata$donor_region=="Multilaterals"], na.rm=T)
})
plot(yrs[yrs>1973], WBcont[yrs>1973], type="l", ylim=c(0, 0.75),
     xlim=c(1974, 2014), xlab="", ylab="% Multilateral Aid from World Bank")


## Figure A2
states <- c("China", "Japan", "Korea South", "Iran", "Spain", 
            "Saudi Arabia", "Mexico", "Turkey", "Italy", "Nigeria",
            "South Africa", "Belgium", "Peru", "Tanzania/Tanganyka", "Uganda",
            "India", "Greece", "Bangladesh", "Kuwait", "Pakistan",
            "Argentina", "Canada", "Sweden", "Germany", "UK", "USA")
VS.ch <- sapply(states, function(x){
  DF$WBvoteshare[DF$year==2014 & DF$state==x] - DF$WBvoteshare[DF$year==1980 & DF$state==x]
})
names(VS.ch)[names(VS.ch)=="Korea South"] <- "South Korea"
names(VS.ch)[names(VS.ch)=="Tanzania/Tanganyka"] <- "Tanzania"
VS.ch <- sort(VS.ch, decreasing=F)

par(mar=c(5.1, 6.1, 4.1, 2.1))
barplot(VS.ch, horiz=T, xaxt="n", yaxt="n", col="lightblue", xlim=c(-.06, .06))
axis(side=1, at=c(-.06, -.04, -.02, 0, .02, .04, .06))
title(xlab="Change in World Bank Vote Power, 1980-2014")
for(i in 1:length(VS.ch)){
  par(las=2)
  mtext(side=2, at=seq(0.75, 31, length.out=length(VS.ch))[i],
        text=names(VS.ch)[i], cex=0.7, srt=90)
}
abline(v=-.06, lty=3, col="lightgray")
abline(v=-.04, lty=3, col="lightgray")
abline(v=-.02, lty=3, col="lightgray")
abline(v=.02, lty=3, col="lightgray")
abline(v=.04, lty=3, col="lightgray")
abline(v=.06, lty=3, col="lightgray")




###############
## IV models ##
###############

## Define functions for clustering SEs
cl <- function(dat,fm, cluster){ 
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }

clusterVCV <- function(data, fm, cluster1, cluster2=NULL) {
  require(sandwich)
  require(lmtest)
  est.fun <- estfun(fm)
  names(fm$model)
  if("as.factor(Year)" %in% names(fm$model)){
    names(fm$model)[names(fm$model)=="as.factor(Year)"] <- "Year"
  }
  inc.obs <- complete.cases(data[,names(fm$model)])
  if("Year" %in% names(fm$model)){
    names(fm$model)[names(fm$model)=="Year"] <- "as.factor(Year)"
  }
  N  <- dim(fm$model)[1]
  NROW <- NROW(est.fun)
  K  <- fm$rank
  cov <- function(cluster) {
    cluster <- factor(cluster)
    u <- apply(est.fun, 2, function(x) tapply(x, cluster, sum))
    meat <- crossprod(u)/N
    M <- length(levels(cluster))
    dfc <- M/(M-1) * (N-1)/(N-K)
    dfc * NROW/N * sandwich(fm, meat=meat)
  }
  cluster1 <- data[inc.obs,cluster1]
  cov1  <- cov(cluster1)
  if(is.null(cluster2)) {
    return(cov1)
  } else {
    cluster2 <- data[inc.obs,cluster2]
    cluster12 <- paste(cluster1,cluster2, sep="")
    cov2   <- cov(cluster2)
    cov12  <- cov(cluster12)
    covMCL <- (cov1 + cov2 - cov12)
    return(covMCL)
  }
}

## exclude Big 4 + India, France to ensure exogeneity of instrument
rm <- which(DF$state %in% c("Russia", "USA", "UK", "China", "India", "France"))
DF2 <- DF[-rm,] 

## Table 1
library(lfe) 
IV1 <- felm(institutional_prolif ~ GDP_lag_log + GDPpc_lag_log + 
              polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
              aid_given_aiddata_lag_log + 
              time_since_prolif + time_since_prolif2 + time_since_prolif3 +
              NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
              region | (VotePowerMisalignment_lag ~ Mikesell_Formula_Change) | state,
            data=DF2)
summary(IV1) # Table 1, col 1

dat <- DF2[-IV1$na.action,]
coef(IV1)[length(coef(IV1))] * sd(dat$VotePowerMisalignment_lag) # effect of 1 sd shift in VPM

## get clustered F-stat 
fs = lm(VotePowerMisalignment_lag ~ Mikesell_Formula_Change + GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
fn = lm(VotePowerMisalignment_lag ~ GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
lmtest::waldtest(fs, fn, vcov = clusterVCV(dat, fs, cluster1="state"))$F[2] # F-test robust to clustering, 22.614
cl(dat, fs, as.character(dat$state)) # first stage results (Appendix Table A3, col 1)


IV2 <- felm(institutional_prolif ~ WB_portsim_st + portsim_noaid + 
              UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log +  
              polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
              aid_given_aiddata_lag_log + 
              time_since_prolif + time_since_prolif2 + time_since_prolif3 +
              NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
              region | (VotePowerMisalignment_lag ~ Mikesell_Formula_Change) | state,
            data=DF2)
summary(IV2) # Table 1, col 2
dat <- DF2[-IV2$na.action,]
coef(IV2)[length(coef(IV2))] * sd(dat$VotePowerMisalignment) # effect of 1 sd shift
coef(IV2)[length(coef(IV2))] * (DF$VotePowerMisalignment[DF$state=="Japan" & DF$year==2012] - DF$VotePowerMisalignment[DF$state=="China" & DF$year==2012]) # effect of giving China Japan's VPM

# f-stat
fs = lm(VotePowerMisalignment_lag ~ Mikesell_Formula_Change + WB_portsim_st + portsim_noaid + 
          UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
fn = lm(VotePowerMisalignment_lag ~ WB_portsim_st + portsim_noaid + 
          UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
lmtest::waldtest(fs, fn, vcov = clusterVCV(dat, fs, cluster1="state"))$F[2] # robust F stat
cl(dat, fs, as.character(dat$state)) # first stage results (Appendix Table A3, col 2)

# create interaction terms for formula shift and state preferences
DF2$Mikesell_Formula_Change_WB_portsim_st <- DF2$Mikesell_Formula_Change * DF2$WB_portsim_st
DF2$Mikesell_Formula_Change_UN_idealpt_USsim_lag <- DF2$Mikesell_Formula_Change * DF2$UN_idealpt_USsim_lag

# first try using these as instruments
IV.i1 <- felm(institutional_prolif ~ WB_portsim_st + portsim_noaid + 
                UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                aid_given_aiddata_lag_log + 
                time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
                region | (VotePowerMisalignment_lag | I(VotePowerMisalignment_lag*WB_portsim_st) | 
                            I(VotePowerMisalignment_lag*UN_idealpt_USsim_lag) ~ Mikesell_Formula_Change + 
                            Mikesell_Formula_Change_WB_portsim_st + Mikesell_Formula_Change_UN_idealpt_USsim_lag) | state,
              data=DF2)
condfstat(IV.i1) # instruments are too weak to yield reliable estimates         

## Instead, use reduced-form regressions

RF <- lm(institutional_prolif ~ Mikesell_Formula_Change + WB_portsim_st + portsim_noaid + 
              UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
              polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
              aid_given_aiddata_lag_log + 
              time_since_prolif + time_since_prolif2 + time_since_prolif3 +
              NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
              region, data = DF2)
dat <- DF2[-RF$na.action,]
cl(dat, RF, as.character(dat$state)) # Table 1, col 3

RF2 <- lm(institutional_prolif ~ Mikesell_Formula_Change + 
               Mikesell_Formula_Change*WB_portsim_st + 
               Mikesell_Formula_Change*portsim_noaid + 
               Mikesell_Formula_Change*UN_idealpt_USsim_lag + 
               WB_portsim_st + portsim_noaid + 
               UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
               polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
               aid_given_aiddata_lag_log + 
               time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
               region,
             data = DF2)
dat <- DF2[-RF2$na.action,]
cl(dat, RF2, as.character(dat$state)) # Table 1, col 4


##########################
## Fixed Effects Models ##
##########################

## Match with CBPS 
library(CBPS)

M <- lm(institutional_prolif ~ VotePowerMisalignment_lag + 
          WB_portsim_st + portsim_noaid + UN_idealpt_USsim_lag + 
          GDP_lag_log + GDPpc_lag_log + 
          region + time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NMC_lag + polity_lag + DBs_joined + 
          aid_received_aiddata_lag_log + aid_given_aiddata_lag_log +
          as.factor(state) + as.factor(year),
         data=DF)
dat <- DF[-M$na.action,]
dat2 <- dat[,c("VotePowerMisalignment_lag", "WB_portsim_st", "portsim_noaid", "UN_idealpt_USsim_lag",
               "GDP_lag_log", "GDPpc_lag_log",
               "NMC_lag", "polity_lag", "aid_received_aiddata_lag_log",
               "aid_given_aiddata_lag_log")]
fit <- CBPS(VotePowerMisalignment_lag ~ WB_portsim_st + portsim_noaid + UN_idealpt_USsim_lag + 
              GDP_lag_log + GDPpc_lag_log + 
              NMC_lag + polity_lag + aid_received_aiddata_lag_log + 
              aid_given_aiddata_lag_log, data=dat2)

## Plot improvement in balance (Appendix Figure A5)
plot(0:1, 0:1, cex=0, xlab="Pearson Correlation",
     yaxt="n", ylab="", xlim=c(-0.7, 0.7))
axis(side=2, at=c(1/4, 3/4), labels=c("CBPS Weighted", "Unweighted"))
abline(v=0, lty=2)
points(balance(fit)[[1]], rep(1/4, 9), pch=19)
points(balance(fit)[[2]], rep(3/4, 9), pch=19)

## Table 2
Match1 <- lm(institutional_prolif ~ VotePowerMisalignment_lag + 
               WB_portsim_st + portsim_noaid + UN_idealpt_USsim_lag + 
               GDP_lag_log + GDPpc_lag_log + 
               region + time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NMC_lag + polity_lag + DBs_joined + 
               aid_received_aiddata_lag_log + aid_given_aiddata_lag_log +
               as.factor(state) + as.factor(year),
          data=dat, weights = fit$weights)
cl(dat, Match1, as.character(dat$state))[c(1:10),] # Table 2, col 1


Match2 <- lm(institutional_prolif ~ VotePowerMisalignment_lag + 
           VotePowerMisalignment_lag*WB_portsim_st + VotePowerMisalignment_lag*UN_idealpt_USsim_lag + 
           WB_portsim_st + portsim_noaid + UN_idealpt_USsim_lag + 
           GDP_lag_log + GDPpc_lag_log + 
           region + time_since_prolif + time_since_prolif2 + time_since_prolif3 +
           NMC_lag + polity_lag + DBs_joined + 
           aid_received_aiddata_lag_log + aid_given_aiddata_lag_log +
           as.factor(state) + as.factor(year),
         data=dat, weights = fit$weights)
cl(dat, Match2, as.character(dat$state))[c(1:10,235:241),] # Table 2, col 2

Match3 <- lm(institutional_prolif ~ VotePowerMisalignment_lag + 
               VotePowerMisalignment_otherDBs_lag + 
               VotePowerMisalignment_lag*WB_portsim_st + VotePowerMisalignment_lag*UN_idealpt_USsim_lag + 
               WB_portsim_st + portsim_noaid + UN_idealpt_USsim_lag + 
               GDP_lag_log + GDPpc_lag_log + 
               region + time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NMC_lag + polity_lag + DBs_joined + 
               aid_received_aiddata_lag_log + aid_given_aiddata_lag_log +
               as.factor(state) + as.factor(year),
             data=dat, weights = fit$weights)
dat2 <- dat[-Match3$na.action,]
cl(dat2, Match3, as.character(dat2$state))[c(1:10,238:242),] # Table 2, col 3



######################
## Robustness Tests ##
######################

## Sensitivity analysis procedure described in "Plausibly Exogenous" by Conley, Hansen, Rossi

delta.range <- seq(from=-3, to=3, by=0.01)
CI.95 <- CI.90 <- as.data.frame(matrix(nrow = length(delta.range), ncol =2))
colnames(CI.95) <- colnames(CI.90) <- c("lower", "upper")
DF2$state <- as.character(DF2$state)
library(ivreg)
library(ivpack)
for(i in 1:length(delta.range)){
  ## re-estimate IV regression with Y-Z(delta) as outcome
  DF2$Sens.Outcome <- DF2$institutional_prolif - (DF2$Mikesell_Formula_Change * delta.range[i])
  model.sens <- ivreg(Sens.Outcome ~ VotePowerMisalignment_lag + 
                        WB_portsim_st + portsim_noaid + 
                        UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                        polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                        aid_given_aiddata_lag_log + 
                        time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                        NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                        region | Mikesell_Formula_Change + 
                        WB_portsim_st + portsim_noaid + 
                        UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                        polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                        aid_given_aiddata_lag_log + 
                        time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                        NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                        region, 
                      data = DF2, x = TRUE)
  CI.95$lower[i] <- coef(model.sens)[2] - cluster.robust.se(model.sens, as.character(DF2$state[-model.sens$na.action]))[2,2] * qnorm(0.975)
  CI.95$upper[i] <- coef(model.sens)[2] + cluster.robust.se(model.sens, as.character(DF2$state[-model.sens$na.action]))[2,2] * qnorm(0.975)
  CI.90$lower[i] <- coef(model.sens)[2] - cluster.robust.se(model.sens, as.character(DF2$state[-model.sens$na.action]))[2,2] * qnorm(0.95)
  CI.90$upper[i] <- coef(model.sens)[2] + cluster.robust.se(model.sens, as.character(DF2$state[-model.sens$na.action]))[2,2] * qnorm(0.95)
  print(i)
}
sens.range <- seq(from=-1.29, to = 1.29, by =0.01) # See if 90% confidence intervals exclude 0 for +-1.29
sens.90 <- c(min(CI.90$lower[delta.range %in% sens.range]), max(CI.90$upper[delta.range %in% sens.range]))
sens.90 # they do.  Effect holds if Z has a direct effect on Y that ranges from -1.29 to 1.29



## Table A4: combine proliferators and joiners

IV1b <- felm(prolif_join ~ GDP_lag_log + GDPpc_lag_log + 
               polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
               aid_given_aiddata_lag_log + 
               time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
               region | (VotePowerMisalignment_lag ~ Mikesell_Formula_Change) | state,
             data=DF2)
summary(IV1b) # Table A4, col 1

IV2b <- felm(prolif_join ~ WB_portsim_st + portsim_noaid + 
               UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
               polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
               aid_given_aiddata_lag_log +
               time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
               region | (VotePowerMisalignment_lag ~ Mikesell_Formula_Change) | state,
             data=DF2)
summary(IV2b) # Table A4, col 2

RFb <- lm(prolif_join ~ Mikesell_Formula_Change + WB_portsim_st + portsim_noaid + 
               UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
               polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
               aid_given_aiddata_lag_log + 
               time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
               region, data = DF2)
dat <- DF2[-RFb$na.action,]
cl(dat, RFb, as.character(dat$state)) # Table A4, col 3

RF2b <- lm(prolif_join ~ Mikesell_Formula_Change + 
                Mikesell_Formula_Change*WB_portsim_st + 
                Mikesell_Formula_Change*portsim_noaid + 
                Mikesell_Formula_Change*UN_idealpt_USsim_lag + 
                WB_portsim_st + portsim_noaid + 
                UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                aid_given_aiddata_lag_log + 
                time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                region,
              data = DF2)
dat <- DF2[-RF2b$na.action,]
cl(dat, RF2b, as.character(dat$state)) # Table A4, col 4


## Table A5: Alternative outcome: bilateral aid

IV1 <- felm(aid_given_aiddata_log ~ GDP_lag_log + GDPpc_lag_log + 
              polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
              time_since_prolif + time_since_prolif2 + time_since_prolif3 +
              NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
              region | (VotePowerMisalignment_lag ~ Mikesell_Formula_Change) | state,
            data=DF2)
dat <- DF2[-IV1$na.action,]
summary(IV1) # Table A5, col 1

IV2 <- felm(aid_given_aiddata_log ~ WB_portsim_st + portsim_noaid + 
              UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
              polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
              time_since_prolif + time_since_prolif2 + time_since_prolif3 +
              NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
              region | (VotePowerMisalignment_lag ~ Mikesell_Formula_Change) | state,
            data=DF2)
summary(IV2) # Table A5, col 2

RF <- lm(aid_given_aiddata_log ~ Mikesell_Formula_Change + WB_portsim_st + portsim_noaid + 
           UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
           polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
           time_since_prolif + time_since_prolif2 + time_since_prolif3 +
           NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
           region, data = DF2)
dat <- DF2[-RF$na.action,]
cl(dat, RF, as.character(dat$state)) # Table A5, col 3

RF2 <- lm(aid_given_aiddata_log ~ Mikesell_Formula_Change + 
            Mikesell_Formula_Change*WB_portsim_st + 
            Mikesell_Formula_Change*portsim_noaid + 
            Mikesell_Formula_Change*UN_idealpt_USsim_lag + 
            WB_portsim_st + portsim_noaid + 
            UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
            polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
            time_since_prolif + time_since_prolif2 + time_since_prolif3 +
            NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
            region,
          data = DF2)
dat <- DF2[-RF2$na.action,]
cl(dat, RF2, as.character(dat$state)) # Table A5, col 4


## Table A6: Alternative independent variable: difference between predicted and observed vote power
## First, predict vote power based on economic fundamentals
EF <- lm(WBvoteshare ~ GDP_lag_log + 
           GDPpc_lag_log + 
           tot.trade + 
           aid_received_aiddata_log + 
           aid_given_aiddata_log, data=DF)
summary(EF)
dat <- DF[-EF$na.action,]
EFp <- predict(EF)
dat$state_year <- paste(dat$state, dat$year, sep="_")
DF$state_year <- paste(DF$state, DF$year, sep="_")

DF$predicted_WBvoteshare <- EFp[match(DF$state_year, dat$state_year)]
cor(DF$predicted_WBvoteshare, DF$WBvoteshare, use="complete.obs")
DF$predicted_WBvoteshare_diff <- DF$WBvoteshare - DF$predicted_WBvoteshare
cor(DF$predicted_WBvoteshare_diff, DF$VotePowerMisalignment, use="complete.obs")

rm <- which(DF$state %in% c("Russia", "USA", "UK", "China", "India", "France"))
DF2 <- DF[-rm,]
IV1c <- felm(prolif_join ~ GDP_lag_log + GDPpc_lag_log + 
               polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
               aid_given_aiddata_lag_log + 
               time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
               region | (predicted_WBvoteshare_diff ~ Mikesell_Formula_Change) | state,
             data=DF2)
summary(IV1c) # Table A6, col 1
dat <- DF2[-IV1c$na.action,]
fs = lm(predicted_WBvoteshare_diff ~ Mikesell_Formula_Change + GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + region + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
fn = lm(predicted_WBvoteshare_diff ~ GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + region + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
lmtest::waldtest(fs, fn, vcov = clusterVCV(dat, fs, cluster1="state"))$F[2] # 5.872227


IV2c <- felm(prolif_join ~ WB_portsim_st + portsim_noaid + 
               UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
               polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
               aid_given_aiddata_lag_log + 
               time_since_prolif + time_since_prolif2 + time_since_prolif3 +
               NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined | 
               region | (predicted_WBvoteshare_diff ~ Mikesell_Formula_Change) | state,
             data=DF2)
summary(IV2c) # Table A6, col 2
dat <- DF2[-IV2c$na.action,]

fs = lm(predicted_WBvoteshare_diff ~ Mikesell_Formula_Change + WB_portsim_st + portsim_noaid + 
          UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + region + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
fn = lm(predicted_WBvoteshare_diff ~ WB_portsim_st + portsim_noaid + 
          UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
          polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
          aid_given_aiddata_lag_log + region + 
          time_since_prolif + time_since_prolif2 + time_since_prolif3 +
          NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined, data = dat)
lmtest::waldtest(fs, fn, vcov = clusterVCV(dat, fs, cluster1="state"))$F[2] # 7.4998



## Table A7: Placebo tests

IV.security <- ivreg(prolif_security ~ VotePowerMisalignment_lag + 
                       WB_portsim_st + portsim_noaid + 
                       UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                       polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                       aid_given_aiddata_lag_log + 
                       time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                       NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                       region | Mikesell_Formula_Change + 
                       WB_portsim_st + portsim_noaid + 
                       UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                       polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                       aid_given_aiddata_lag_log +
                       time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                       NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                       region, 
                     data = DF2)
dat <- DF2[-IV.security$na.action,]
cl(dat, IV.security, as.character(dat$state)) # Table A7, col 1

IV.trade <- ivreg(prolif_trade ~ VotePowerMisalignment_lag + 
                    WB_portsim_st + portsim_noaid + 
                    UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                    polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                    aid_given_aiddata_lag_log + 
                    time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                    NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                    region | Mikesell_Formula_Change + 
                    WB_portsim_st + portsim_noaid + 
                    UN_idealpt_USsim_lag + GDP_lag_log + GDPpc_lag_log + 
                    polity_lag + NMC_lag + DBs_joined + aid_received_aiddata_lag_log + 
                    aid_given_aiddata_lag_log +
                    time_since_prolif + time_since_prolif2 + time_since_prolif3 +
                    NatIncome_Combined + Imports_Combined + Exports_Combined + Gold_FEX_Combined + 
                    region, 
                     data = DF2)
dat <- DF2[-IV.trade$na.action,]
cl(dat, IV.trade, as.character(dat$state)) # Table A7, col 2


## Figure A3: Covariate Balance for States affected by Mikesell Formula
## Create balance plot
DF$VoteShift_gain <- ifelse(DF$Mikesell_Formula_Change > 0, 1, 0)
DF.sub <- DF[-which(DF$state %in% c("China", "Russia", "USA", "UK", "France", "India")),]
covs <- DF.sub[DF.sub$year == 1947,colnames(DF.sub) %in% c("GDP_lag_log", "NMC_lag", "polity_lag", 
                                                           "GDPpc_lag_log", "aid_received_aiddata_lag_log", "year",
                                                           "aid_given_aiddata_lag_log", "WB_voteshare", "VoteShift_gain",
                                                           "WBvoteshare")]
GDP.diff <- mean(covs$GDP_lag_log[covs$VoteShift_gain==1], na.rm=T) - 
  mean(covs$GDP_lag_log[covs$VoteShift_gain==0], na.rm=T)
GDP.diff.st <- GDP.diff / sd(covs$GDP_lag_log, na.rm=T)
GDPt <- t.test(covs$GDP_lag_log[covs$VoteShift_gain==1] / sd(covs$GDP_lag_log, na.rm=T),
       covs$GDP_lag_log[covs$VoteShift_gain==0] / sd(covs$GDP_lag_log, na.rm=T))
GDPpc.diff <- mean(covs$GDPpc_lag_log[covs$VoteShift_gain==1], na.rm=T) - 
  mean(covs$GDPpc_lag_log[covs$VoteShift_gain==0], na.rm=T)
GDPpc.diff.st <- GDPpc.diff / sd(covs$GDPpc_lag_log, na.rm=T)
GDPpct <- t.test(covs$GDPpc_lag_log[covs$VoteShift_gain==1] / sd(covs$GDPpc_lag_log, na.rm=T),
       covs$GDPpc_lag_log[covs$VoteShift_gain==0] / sd(covs$GDPpc_lag_log, na.rm=T))
polity.diff <- mean(covs$polity_lag[covs$VoteShift_gain==1], na.rm=T) - 
  mean(covs$polity_lag[covs$VoteShift_gain==0], na.rm=T)
polity.diff.st <- polity.diff / sd(covs$polity_lag, na.rm=T)
polityt <- t.test(covs$polity_lag[covs$VoteShift_gain==1] / sd(covs$polity_lag, na.rm=T),
       covs$polity_lag[covs$VoteShift_gain==0] / sd(covs$polity_lag, na.rm=T))
NMC.diff <- mean(log(covs$NMC_lag[covs$VoteShift_gain==1]), na.rm=T) - 
  mean(log(covs$NMC_lag[covs$VoteShift_gain==0]), na.rm=T)
NMC.diff.st <- NMC.diff / sd(log(covs$NMC_lag), na.rm=T)
NMCt <- t.test(log(covs$NMC_lag[covs$VoteShift_gain==1]) / sd(log(covs$NMC_lag), na.rm=T),
       log(covs$NMC_lag[covs$VoteShift_gain==0]) / sd(log(covs$NMC_lag), na.rm=T))
WBv.diff <- mean(covs$WBvoteshare[covs$VoteShift_gain==1], na.rm=T) - 
  mean(covs$WBvoteshare[covs$VoteShift_gain==0], na.rm=T)
WBvt <- t.test(covs$WBvoteshare[covs$VoteShift_gain==1],
       covs$WBvoteshare[covs$VoteShift_gain==0])

par(mar=c(5.1,8,4.1,2.1))
plot(c(GDP.diff.st, GDPpc.diff.st, NMC.diff.st, polity.diff.st), 
     c(2,1,0,-1), xlim = c(-2.5, 3), pch=19,
     yaxt = "n", xlab = "", bty="n", xaxt="n",
     ylab = "", ylim = c(-1.1, 2))
abline(v=0, lty=2)
lines(GDPt$conf.int, c(2,2))
lines(GDPpct$conf.int, c(1,1))
lines(NMCt$conf.int, c(0,0))
lines(polityt$conf.int, c(-1,-1))
axis(side=2, at = -1:2, las=1, pos=-2,
     labels = c("Polity", "Military Power", 
                "GDP per capita", "GDP"))
axis(side=1, at = seq(from=-2, to=2, by=0.5))
text(2.05, -1, paste("p =", as.numeric(round(polityt$p.value, 2))))
text(2.05, 0, paste("p =", as.numeric(round(NMCt$p.value, 2))))
text(2.05, 1, paste("p =", as.numeric(round(GDPpct$p.value, 2))))
text(2.05, 2, paste("p =", as.numeric(round(GDPt$p.value, 2))))






